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3-d resistive MHD simulations of magnetic reconnection and the tearing mode instability in 

current sheets 
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Magnetic reconnection plays a critical role in many astrophysical processes where high energy emission is observed, 
e.g. particle acceleration, relativistic accretion powered outflows, pulsar winds and probably in dissipation of Poynting 
flux in GRBs. The magnetic field acts as a reservoir of energy and can dissipate its energy to thermal and kinetic 
energy via the tearing mode instability. We have performed 3d nonlinear MHD simulations of the tearing mode 
instability in a current sheet. Results from a temporal stability analysis in both the linear regime and weakly nonlinear 
(Rutherford) regime are compared to the numerical simulations. We observe magnetic island formation, island merging 
and oscillation once the instability has saturated. The growth in the linear regime is exponential in agreement with 
linear theory. In the second, Rutherford regime the island width grows linearly with time. We find that thermal energy 
produced in the current sheet strongly dominates the kinetic energy. Finally preliminary analysis indicates a P(k) 
4.8 power law for the power spectral density which suggests that the tearing mode vortices play a role in setting up 
an energy cascade. 
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1. Introduction 

Magnetic reconnection plays a critical role in many 
astrophysical processes, e.g. particle acceleration [1], 
accretion disks [2] and solar flares [3] . It is also impor- 
tant in laboratory fusion devices such as tokamaks. 
Magnetic reconnection is a topological change in the 
field which violates the frozen-flux condition of ideal 
magnetohydrodynamics (MHD). If a magnetic field 
can leak across the plasma it can reach a lower energy 
state - in the case of a current sheet it can undergo 
"tearing" into filaments or magnetic islands. A cur- 
rent layer of thickness a may dissipate on timescales 
shorter than the resistive timescale a 2 /rj due to the 
tearing mode instability (hereafter TMI). The TMI 
was first discovered in tokamaks and stellerators and 
has been extensively studied since the pioneering 
works [4] and [5] (see also [6], [7], [8]). In this paper 
we present numerical simulations of the dissipation 
of current sheets and the formation of magnetic is- 
lands due to the tearing mode instability. In Section 
2 we remind the reader of the predictions of the linear 



and weakly nonlinear analyses of Furth et al. [4] and 
Rutherford [5]. In Section 3 we present our numeri- 
cal method. In Section 4 we compare the results with 
linear and weakly nonlinear theory. In Section 5 we 
present power spectra derived from 3d simulations. 
In Section 6 we discuss the pitfalls of simulations 
where reconnection is not properly tracked and the 
implications for reconnection-driven turbulence. 



2. Theory 

2.1. The linear regime 

In the linear analysis it is assumed that resistivity 
is only important within the current layer. Separate 
solutions may then be derived for the exterior, ideal 
MHD region and the interior, resistive region. Us- 
ing asymptotic matching [4] derived stability range 
and growth rates for the TMI in the inviscid linear 
regime, neglecting the effects of compressibility. The 
range of the instability is a < 1, a = ak, where k is 
the wave number and a is the current sheet width. 
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The growth rate is: 
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(1) 



The linear growth rate is a function of the mode 
number, k, the current sheet width, a, and the 
Lundquist or magnetic Reynolds number, S. The 
equations are valid for a > S'~*. The system stays 
linear for as long as the magnetic island does not 
exceed the width of the current sheet. 

2.2. Non-linear regime 

Rutherford [5] studied the evolution of the TMI in 
the weakly nonlinear regime, by considering the ef- 
fect of second order eddy currents on the current 
sheet. [5] found in the nonlinear regime that the is- 
land width grows as t, instability growth rate slows 
down from exponential to t 2 , and that the critical 
amplitude where the linear solution ceases to be valid 
is: 



a Harris current sheet of the form B y = tanh(ir). 
We initialise the field in the magnetic vector po- 
tential A z = ln|cosh(x)| and take the curl to de- 
rive B. To this we add a perturbation of the form 
cos(kx) in 2d and cos(kx)cos(ky) in 3d. Free param- 
eters are a, (3 = -J^ and the resistivity r\. Zero- 
gradient Neumann-like "outflow" boundary condi- 
tions often produce subsonic reflection for compress- 
ible flows [10, 11]. We use a coarse mesh boundary 
which reduces this effect. 

4. Simulation Results 
4.1. 2D 



(2) 





Table 1. Growth Rates and Saturation Times for 
Lundquist number 2400 
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Fig. 1. Formation of magnetic islands: magnetic field lines 
and density colourmap are shown at times: 245, 365, 490. 




2.3. Previous work 

More recent work has concentrated on including the 
physics of the Hall effect using multifluid or Hall 
MHD and electron inertia, using PIC techniques. 
Significant increase in reconnection rates is found in 
these studies [9] . However in these works the authors 
do not compare against analytical theory nor is the 
3D power spectrum calculated. 

3. Numerical setup 

We perform first 2D and then 3D simulations in re- 
sistive MHD using the astrophysical code PLUTO. 

As in any numerical code some numerical resis- 
tivity is present. In the linear analysis numerical re- 
sistivity has been assumed to be small in comparison 
with the physical resistivity. Our initial mean field is 



Fig. 2. Same as Figure 1 but for a = 0.3 

We performed 2d and 3d simulations of cur- 
rent sheets for different initial perturbations k = 
0.01, 0.04, 0.3, 0.6. We tracked the growth rate of the 
cross-sheet magnetic field B x . In the simulations an 
initial transient period, t tra nsient ~ lOOr^ roughly 
corresponding to half the crossing time of the do- 
main {t = 112ta) was observed before the linear 
mode was established. As the tearing mode grows 
linearly, field lines reconnect, new separatrices are 
formed and magnetic islands grow in size, eventually 
exceeding the size of the original current sheet. In 
Figures 1 and 2 the time evolution of magnetic is- 
lands for a = 0.01,0.3 are plotted against time in 
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units of the Alfven crossing time, ta- 

2d Growth of log e lB x mox l cc = 0.01 




2000 



Fig. 3. Growth of log of maximum B x plotted against time. 
The linear growth rate is plotted and the critical value for 
saturation. 



2d Growth of log e lB x mox l a= 0.3 
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Fig. 4. Same as Fig. 3 but for a = 0.3. 



ear regime is found with slopes near the analytical 
predictions of Furth et al. [4]. 

Island width evolution a = 0.01 




500 



Fig. 6. Magnetic island width evolution plotted against time 
in units of ta for a = 0.01. The width is estimated using the 
full-width half maximum of the thermal pressure. It remains 
approximately constant in the linear phase and transfers to a 
linear growth in t once in the Rutherford regime. 



Island width evolution cx = 0.30 
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2d Energy Evolution a = 0.01 
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Fig. 5. Time evolution of change in kinetic energy, magnetic 
energy and internal energy. Most of the magnetic energy is 
converted into thermal energy. 



4.1.1. Agreement with linear analysis 

In Figures 3 and 4 the growth rate in cross-sheet 
magnetic field, B x for different exciting modes of 
the TMI is plotted against time. In all cases a lin- 



Fig. 7. Same as 6 but for a = 0.3. 

4.1.2. Agreement with nonlinear analysis 

We find plotting the island width against time that 
it remains small in the linear regime and grows as t 2 
in the Rutherford regime. We can estimate the time 
of transition from linear to Rutherford from both the 
B x plots and the island width plots, however the is- 
land width plots provide a more accurate estimate. 
The time at which the system becomes nonlinear is 
tabulated in Table 1. The magnetic island width can 
be derived from eqn 18. in [5, 6] and can be written 
in approximate form as 



Wi 



t-t d 



(3) 



where td accounts for the total delay (transient and 
time spent in linear regime) in our simulations before 
entering the Rutherford regime. 
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5. 3D Power spectra 

In Figure 8 we plot the power spectrum for (3 — 2 
3d current sheet. The slope of P{k) is —4.8, within 
the range of values found by [12] for /3 = 1 and 
/? = oo. [12] note that in their driven, supersonic 
turbulence simulation they have a constant beta and 
a constant mass-to-flux ratio "modulo the numeri- 
cal reconnection effect". We obtain the same power 
spectrum scaling without any added velocity pertur- 
bations. 
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Fig. 8. 3d kinetic power spectrum at time t = 295ta with 
ovcrplottcd fc~ 4 ' 8 

6. Discussion and Conclusion 

6.1. Consequences for driv en-turbulence 
simulations 

Our results show that even without either a driven 
turbulence mechanism, or an initially turbulent ve- 
locity spectrum it is possible for the vortices gener- 
ated by the TMI to mimic the power spectrum seen 
in simulations of turbulence. This can have serious 
consequences for MHD simulations where the resis- 
tivity is not explicitly constrained (i.e. numerical re- 
sistivity plays a role), since there will be inevitably be 
some numerical-reconnection driven vortices present 
in the simulation. However using an explicit resistiv- 
ity, the contribution from reconnection can be esti- 
mated using the formulae in [4, 5]. 

Finally, we propose as a benchmark for resis- 
tive MHD codes the tearing mode instability test, 
which can be compared with both analytical results 
for both the linear and nonlinear regimes, as well 



as with laboratory experiments. The test may also 
prove useful as a way of quantifying the effects of 
numerical resistivity in an MHD code. 
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